A Full-Matrix Approach for Solving General Plasma Dispersion Relation 
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A hitherto difficult and unsolved issue in plasma physics is how to give a general numerical solver 
for complicated plasma dispersion relation, although we have long known the general analytical 
forms. We transform the task to a full-matrix eigenvalue problem, which allows to numerically 
calculate all the dispersion relation solutions exactly free from convergence problem and give po- 
larizations naturally for arbitrarily complicated multi-scale fluid plasma with arbitrary number of 
components. Attempt to kinetic plasma via A''-point Pade approximation of plasma dispersion 
function also shows good results. 



Since only few simple dispersion relations are analyt- 
ical tractable in plasma physics, it is a historical issue 
to develop general numerical solvers for practical ap- 
plications, especially for space, astrophysical and laser 
plasma. However, it is a hitherto difficulty to develop 
a kinetic solver with general distribution function and 
general other effects. The first problem comes from the 
Landau contour integral (e.g.. Landau damping) of the 
kinetic distribution function. Second problem is the in- 
finity orders of Bessel function summation for magnetized 
plasma which is related to the cyclotron resonance (e.g., 
Bernstein modes). These two obstacles make it difficult 
to develop a root finding solver with good convergence. 

WHAMP (Waves in Homogeneous Anisotropic Multi- 
component Magnetized Plasma) code by Ronnmark[Tl[5] 
is an important step to that goal, which can be used 
to calculate general non-relativistic kinetic wave disper- 
sion relation in plasmas with parallel beams. To calcu- 
late fast, Pade approximation is used for plasma disper- 
sion function. However, the numerical convergence is still 
not as expectations, especially that it is difficult to give 
proper initial guesses for high frequency (e.g., ui > lOftci) 
modes. 

Bret [3] discussed how to derive and solve the 3-by-3 
parallel beam-plasma dielectric tensor in fluid approxi- 
mation with relativistic effect with the help of computer 
(Mathematica). However, it is still not easy to use this 
method if there are many components. 

In our treatment, we do not need to derive the fi- 
nal 3-by-3 dispersion relation tensor matrix for E = 
(Ex, Ey, Ez) as the conventional treatment, such as by 
Stix[6^ and used by Ronnmark[Tl [2] and Bret et oL[3H5], 
which is helpful to provide analytical insight but not a 
must for numerical solver. 

In this manuscript (brief communication), as the first 
topic, we present a trivial multi-fluid non-relativistic 
magnetized arbitrary orient warm beam plasma problem, 
which is very difficult via usual treatments, to show how 
our full-matrix approach works. 
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The original equations are 

dtUs = -V • (nsVs), 

dtE ^c^V X B- J/eo, 
dtB = -V xE, 



(la) 
(lb) 
(Ic) 
(Id) 



with 



J ^J2s1>^'^sVs, (2a) 

dtiPsPj''^) ^ 0, (2b) 

where ps = rrisUs, Ps = kBTso/ms and c^ = 1/p.otQ. 

In cold plasma limit, for parallel beam, the usual 3- 
by-3 dispersion relation tensor D{k,uj) can be derived 



as 
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and n = ck/uj, lo'^ = u — k ■ Vso, Bq = (0, 0, i?o), 
VsQ = (0,0, Uds), k = {kx,0,kz) = {ksm0,O,kcos9), 
^cs = qsBo/ms (note: Qe = -e) and w^^ = n^Qq^J earns. 
The dispersion relations is 



r^-i>. 



D{k,uj) = det[D{k,Lo)]^Q. (5) 

The summation ^^ in Kij [i,j = x, y, z) with w in the 



denominator, which is singularity at w = and uj 



'2 



w^s, makes a general good convergence numerical solver 
very difficult. And, we need also give good initial guesses 
when using special root finding solver such as Newton's 
iterative method or give a guess domain in complex plane 
using such as Davies' method[7J. 

We can get an explicit polynomial form dispersion re- 
lation equation from (Is]) for k{uj,6) easily. While, it is 
very cumbersome to calculate an explicit form for a;(fc, 6), 
even though with the aids of computer (e.g., using Math- 
ematica). Without beam (vds = 0) and for only one ion 
species (s — e, i), a fifth order explicit form polynomial 
for a;^(fc,6') is given in Swanson's textbookjS], which is 
simplified from hundreds of terms. 

The above investigations imply that it is not a satis- 
factory choice to solve the dispersion relation using (Is]) 
directly as usual treatment. 

An alterant method is using the original full dispersion 
relation matrix and then treating the task as a matrix 



eigenvalue problem, which need neither derive the final 
dispersion relation equation nor worry about how to solve 
it. 

The linearized version of ( 11 ) with / — /o + /ie''""^~"^* 



is equivalent to a matrix eigenvalue problem (similar 
treatment can be found in [S] for MHD equations) 



XX ^M-X. 



(6) 



with X — ~iuj the eigenvalue and corresponding X eigcn 
vector, which represents the polarization information of 
each normal/eigen mode solution. 



For (111, X is 



{{nsl,Vslx,Vsly, Vslz}, Eix,Eiy,Eiz, Bix,Biy, Biz) 

and M is 
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TABLE I: Comparing the cold plasma solutions using matrix 
method and Swanson's polynomial. Only positive solutions 
are shown, since oji = uPl. 



c." 


C.S 


UJ^' 


a;S 


10.5152 


10.5152 


1.1330E-4 


1.1330E-4+ilE-16 


10.0031 


10.0031 


- 


1E-32-M1E-18 


9.5158 


9.5158 


- 





2.4020E-4 


2.4020E-4-i3E-17 







where we have used psx = ^sPsoPsil Psa, which is from 
(2b), and c^ = jsPso/ Pso- The effects from non-zero equi- 



librium quantities (Jq and v^o x Bq) are omitted here 
(Note: This is another annoying unsolved problem in lit- 
eratures. Principally, we need treat it as inhomogeneous 
problem.). 

For s kinds of species, the dimensions of M are (4s -I- 
6) X (4s -|- 6). Here, we can get all the solutions of the 
above system exactly (without convergence problem) via 
standard matrix eigenvalue solver, e.g. function eig() in 
MATLAB. 



Cold limit (p^o — 0), without beam {vso — 0),s — e,i, 
numerical solutions of ^ (w^) and Swanson's polyno- 
mial (oj^) are given in Table HJ with kc = 0.1, 9 = 7r/3, 
milrrie = 1836 and Wpg = lOwce. The results are consis- 
tent with each other exactly, except some small (< 10^^^) 
numerical errors. 

A 3-component result with perpendicular electron 
beams is shown in Figfll 

Using this method for other arbitrarily complicated 
fluid system is just straightforward, due to that it is just a 
directly rewriting of original linearized equations without 
any approximations. 

As a second topic, we try this matrix method for ki- 
netic plasma. 

A possible general method is discreting[Tni [H] the 
distribution function in w-space or using basis function 
expansionjllj to transform the equations to matrix form. 
However, it is shown that for both Vlasov-Possion[TD] 
and Vlasov- Ampere [H] systems, the (Landau) damping 
(5(a;) < 0) normal modes are not eigenmodes but just 
related to spectral density accumulating of eigenmodes, 
though those methods can give some correct solutions 
for growth modes. Another problem is that the matrix 
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FIG. 1: Test run for fluid beam plasma, with Cie,2e,i = 0.01c, 
viex = O.Olc, nie = 0.9, V2ex = — 0.09c and n2e = 0.1. 



dimension should be very large to calculate accurately 
since discrete of v is introduced. 

Another attempt is the initial value version[T21 [13] of 
this matrix method, which can give correct kinetic Lan- 
dau damping, but gives only few largest imaginary part 
modes and can not overcome multi-scale problem. 

Here, as our first successful attempt to kinetic plasma, 
we limit to simple multi-component ID electrostatic 
(ESID) problem with drift Maxwellian distribution, with 
dispersion relation 
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TABLE II: Comparing the Landau damping solutions using 
matrix method and original Z{C,) function. 
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FIG. 2: Comparing of all solutions from matrix method and 
Z{() function 



where A|,, - ^ 



nsQs 



Vts 



_ l2kBT, 



and C,s — 



^^^f^^^,and 



use A'^-point Pade approximation of plasma dispersion 
function as by Ronnmark [U |2] 
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where A = 8 is used, which shows to be very accurate [TTl 



for most domain in upper plane (except y < ^fnx^e^^ , 
x 3> 1, with s — X + iy). Bad performances are just for 
strong damping domain, for which we have little interest. 
Combining (Is]) and ^ , gives a very similar dispersion 
relation equation 



EE 
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{UJ ~ CsjY 



0, 



(10) 



as the one of the fluid ESID beam plasma, with hgj 

2\'^" and Csj 

form the problem to an equivalent linear system 



k{vsQ + VtsCj), which can help us trans- 



ujA 



T-sj" — OgjUgj -\~ CgjAgj, 


(11a) 


cjB.j = Cs,B,j + C, 


(lib) 


C^Es.As,, 
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which is an eigenvalue prolem of a 2sN x 2sN dimen- 
sions eigen matrix M, with sN = sxN. The sigularity of 
dominator, which meet in conventional treatment, is can- 
celed after this transformation. And the matrix method 
can support multi-component very easily and naturally. 
For Langmuir wave Landau damping, calculating the 
largest imaginary part solution using matrix method 
(w*^) and original Z{(^) function (w^) are shown in Table 



Fig. [5] shows all the solutions of matrix method and 
also solutions using Z{C,) function for k = 0.8. We can 
see that only the largest imaginary part solutions are 
overlapping, which is enough for what we want. 

For two-frequency-scale ion acoustic mode, besides the 
Langmuir mode to — 2.0458 — 0.8513i, the largest imag- 
inary part solution in matrix method is also consistent 
with Z{C.) function solution, e.g., Ti = Tg, mi = 18367Tie, 
kXoe = 1, gives a; = 0.0420 - 0.0269i. 

First four largest imaginary part solutions of electron 
bump-on-tail mode, with ion effects, are shown in FigIS] 
for Tb ^ Te = Ti, Vb = 5vte and rif, = 0.1. For this 
test run, matrix method calculates very fast and auto- 
matically, and can give all the solutions we want. While, 
using Z{C,) function, we need test different initial guess 
one by one to select different modes (not shown here). 

Detail features and benchmarks of highly nontriv- 



(a) growth rate, n =0.1 (b) frequency, T =T , v =5v 
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FIG. 3: ESID electron bump-on-tail modes with ion effect. 



ial solvers for multi-component magnetized and un- 
magnetized fluid plasma XBX = AX with and without 



anisotropic, relativistic, beam and gradient effects, will 
be reported elsewhere, where B j^ I is mainly from rel- 
ativistic factor. Though haven't shown yet, we hope, in 
the future, this method can also performances well for 
other nontrivial kinetic problems. 

In summary, a very efficacious method is presented to 
solve both fluid and (simple) kinetic plasma dispersion 
relation, which have overcome many (almost all) troubles 
in conventional treatments. 

[The author would like to thank M. Y. Yu for improv- 
ing the manuscript.] Conversations with Liu Chen, Y. 
Xiao and other researchers at IFTS-ZJU are appreciated. 
Information from Y. Lin, X. Y. Wang, R. Denton and 
Ling Chen are also acknowledged. 

Codes for solving ^ and ( 11 ) are provided as supple- 
mentary materials for one who hopes to quickly access to 
this matrix approach. 
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